Phase separation of trapped spin-imbalanced Fermi gases in one-dimensional optical lattices 
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We calculate the density profiles of a trapped spin-imbalanced Fermi gas with attractive interactions in a 
one-dimensional optical lattice, using both the local density approximation (LDA) and density matrix renor- 
malization group (DMRG) simulations. Based on the exact equation of state obtained by Bethe ansatz, LDA 
predicts that the gas phase separates into shells with a partially polarized core and fully paired wings, where the 
latter occurs below a critical spin polarization. This behavior is also seen in numerically exact DMRG calcula- 
tions at sufficiently large particle numbers. Unlike the continuum case, we show that the critical polarization is 
a non monotonic function of the interaction strength and vanishes in the limit of large interactions. 
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I. INTRODUCTION 

The physics of population imbalanced Fermi gases in one 
dimension is currently attracting substantial interest since in 
these systems, the one-dimensional (ID) counterpart of the 
Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state [lLl2l] can be 
realized. This inhomogeneous superfluid state of fermions 
was suggested long time ago, yet its unambiguous observa- 
tion in condensed matter systems has turned out to be quite 
challenging lUrB]. Ultracold atomic gases offer a perfect av- 
enue to search for exotic pairing states such as the FFLO state, 
due to the clean experimental conditions and the tunability 
of interactions and imbalance ■ First experiments for 
three-dimensional (3D) systems at by Ketterle and co-workers 
ifTolfTlll and by Hulet and co-workers lfl2j[T3ll as well as theo- 
retical work for 3D lfl4ll indicate that the chances of observing 
the FFLO state in 3D are slim since this phase appears to be 
stable only in a narrow window of the BCS-BEC crossover 
phase diagram. 

A promising direction to detect the FFLO state is to confine 
atoms in ID configurations. In this limit the Fermi surface 
reduces to two points, kpf and kpi, forcing all Cooper pairs 
to share the same wave-vector Qfflo = &F,t — ^f,i- For 
systems of fermions with attractive interactions - usually de- 
scribed by the Hubbard model or the associated continuum 
model - it has indeed been shown by means of analytical and 
numerical techniques that any finite imbalance drives the sys- 
tem into the ID FFLO state, which has algebraically decaying 
pair-pair correlations in the s-wave channel, modulated with 
the FFLO wavevector Qfflo (HUH- 

In the more general case of a two-channel model that also 
accounts for the formation of composite molecules, the FFLO 
phase still emerges on the BCS side in the imbalanced case, 
yet eventually gives room to a Bose-Fermi mixture phase 
l26ll . The FFLO phase also exists in mass-imbalanced 



systems in ID (2 
tures (see Ref. ||32 



■ 31|], which could be realized in Li-K mix- 
]). It should further be possible to drive 
the FFLO state in a globally balanced mixture by using spin- 
dependent potentials yH or spin-dependent optical lattices 
I134tl . Furthermore, FFLO physics has been studied in 3D ar- 
rays of coupled ID systems ll23l l24l r35ll . the two-leg ladder 



geometry 13611 . and the possibility of inhomogeneous pairing 
states in multi-band systems has been explored in Ref. l37ll . 
Finally, we mention that the spectral function of attractively 
interacting fermions with spin-imbalance has been calculated 
and discussed in Ref. 0811 . 

Experimentally, the simplest observables that can be mea- 
sured in trapped Fermi gases are the density profiles ifiol - 
[HH1] of the two spin components iflO^ll [39ll . On the basis 
of the exact Bethe ansatz solution for the homogeneous sys- 
tem, it has been independently predicted by Hu et al. lfl6ll and 
one of the authors 111711 that the trapped gas phase-separates 
into two shells: a partially polarized core, and either fully 
paired or fully polarized wings. Very recently, this theoretical 
prediction has been addressed and verified in an experiment 
with ultracold fermions in ID wave-guides at Rice 114011 - 

In the case of an optical lattice, the shell structure of spin- 
imbalanced fermions was studied with the density matrix 
renormalization group (DMRG) method in Refs. lfl8l[T9l liHl 
and with quantum Monte Carlo (QMC) in Ref. l20ll. In 
Ref. lfl8ll by some of us, clear evidence for the presence of 
FFLO correlations of the partially polarized core has been re- 
ported on systems with up to 80 fermions. In another DMRG 
work by Tezuka and Ueda lfl9ll . a phase diagram for the shell 
structure has been proposed, with lower critical polarizations 
for the emergence of the fully polarized wings. In Ref. Ql, 
the presence of the fully polarized wings was deduced from 
the presence of certain local maxima in the density difference 
profile in the outer regions of the trap. However, neither of 
these studies observed the fully paired wings, which have been 
seen in QMC calculations, in the continuum case l22tl . 

The purpose of this work is to clarify this open point, 
namely the presence of fully paired wings in an optical lat- 
tice. As a result, numerically exact DMRG simulations indeed 
show the existence of fully paired wings at low spin polariza- 
tions. To that end we perform large scale DMRG l43l liill 
calculations of the ID Hubbard model with (effectively) up to 
320 fermions in a harmonic trap, which is comparable to the 
particle numbers typically encountered in experiments lp7L l45fl . 
We carefully compare the DMRG results with those obtained 
starting from the exact solution of the homogeneous lattice 
gas, including the effect of the trap via the local density ap- 
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FIG. 1: (color online) Grand-canonical phase diagram of the ID 
attractive Hubbard model l[T) for U = — 5f (compare (42tl ). The 
vertical lines are LDA trajectories yielding different shell structures 
for the trapped gas: the dashed (red) and the solid (blue) lines cor- 
respond to configurations with fully paired (ED) and fully polarized 
(FPi) wings, respectively. The dotted line corresponds to a special 
value of the magnetic field, h = \/U 2 + 16t 2 /2 — 2t, where the 
partially polarized phase extends to the edge of the cloud. 



proximation (LDA), as previously done for the continuum 
case QSm. 

Our study unveils that the fully paired wings are present be- 
low a critical spin polarization, which we calculate as a func- 
tion of effective density and interaction strength. In contrast 
with the continuum case, in the strong coupling regime, the 
fully paired wings shrink in size as the attraction strength in- 
creases. As we shall see, this effect is due to the enhanced 
inertia of the molecules with respect to unpaired fermions. 

The outline of this work is as follows. In Sec. [Til we shall 
introduce the model and discuss in detail how we obtain the 
density profiles of the trapped gas, using either DMRG sim- 
ulations or LDA applied to the exact equation of state for 
the homogeneous system. Our main results are contained in 
Sec. [Hi] where we compare the two methods and discuss den- 
sity profiles as well as the critical polarization. 



II. MODEL, METHODS AND SET-UP 
A. Hamiltonian 

Our simulations are performed for the ID Hubbard model 
with an attractive onsite interaction U < 0: 

L-l L 

H = -t ^ [c\ a Ci+ia + h.c)j + U n^nn (1) 

i— 1,(7 i—1 

where c\ a creates a fermion with spin a =t,| at site /, 
?i£o- = c^cta, rig = n,£-f + ri(i is the local density, and t 
is the hopping matrix element. We define x = ia, where a 
is the lattice spacing, set to unity. The local polarization is 
defined as (s;) = (n,-^ — n,j,), the filling factor is n = N/L 



where N a = an d N — N a . Analogously, we 

define the densities for spin up and down as n a = N a /L, 

Our main interest will be in the case of trapped fermions for 
which we can write the Hamiltonian as 

H = Hq + i?trap (2) 
L 

tftrap = vY / (^-L/2) 2 n t , (3) 

i=l 

where in the second line, V is the trapping potential. To char- 
acterize the trapped system, a suitable choice of variable is the 
total polarization p = (Nf — Ni)/N and the effective density 
Pcft = NW. In the following we fix the energy scale by 
setting t = 1. 



B. Bethe ansatz and Local Density Approximation 

The density profiles of the two component Fermi gas can be 
obtained starting from the exact Bethe ansatz solution of the 
homogeneous system ([T), taking into account the effect of the 
trap via the local density approximation. 



1. Homogeneous system 

The exact ground state energy per site E(n 
in the thermodynamic limit is given by (|4; 
erences therein) 
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where pi(k) and p2(k) are spectral functions satisfying the 
following coupled integral equations: 
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(5) 



Here K v (x) = 2vu/(v 2 u 2 + x 2 ), with u = \U\/A. The inte- 
gration limits B, Q are related to the total density n = n-f + n\. 
and the density difference s = n-j- — by the conditions 



pi(k)dk = s, 



p 2 {X)dX = n — s. (6) 



Equations (lUi-© are only valid below unit filling, n < 1. 
For higher concentrations the ground state energy can be ob- 
tained through a particle-hole transformation h a = and is 
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FIG. 2: (color online) Particle density (rii) (solid lines) and local 
polarization (si) (dashed lines) for U = -8t, V = 0.0005i, and (i) 
a symmetric trap with N = 80 (lines) and (ii) an asymmetric trap 
with N = 40 (lines with symbols), (a) p = 0; (b) p = 1/2. 



given by s) = £^(2 — n, s) + U(n — 1). From the en- 
ergies, we calculate the averaged chemical potential and the 
effective magnetic field: 



dE_ 

dn ' 



dE 
ds 



(1) 



Interpreting Eq. (|7]i as a change of variables, one obtains the 
phase diagram of the Hamiltonian (Q~|), as shown in FigJT|for a 
specific value of the interaction U = —5t (compare 1H2I0 . We 
can distinguish different phases for a low filling n = N/L < 
1: (i) the vacuum (V), corresponding to s = n = 0, (ii) the 
unpolarized phase (ED) characterized by equal densities of 
majority and minority fermions, (iii) the partially polarized 
phase (PP), corresponding to < s < max(n, 2 — n); (iv) 
the fully polarized phase with n < 1 (FPi), and (v) the fully 
polarized phase with n = 1 (FP2). The remaining phases at 
higher densities correspond to situations where at least one of 
the component forms a band insulator. 



2. Trapped system 

The density profiles of the gas in a shallow trap can be cal- 
culated via the local density approximation, assuming that the 
system is locally homogeneous: 



fj,[n(x), s(x)] 
h[n(x), s(x)] 



h°, 



Vx z 



(8) 
(9) 



where /i° and h° are the chemical potential and the magnetic 
field at the trap center. In the second line of Eq. ©, we have 
exploited that the external potential is the same for both com- 
ponents, corresponding to vertical trajectories in the phase di- 
agram of Fig.Q] The constants fi°, h° are fixed by the overall 
spin populations N a through the normalization conditions 



(10) 



N = n(x)dx, N t - N±= / s(x)dx 



By inverting Eq. (0 we see that the normalization conditions 
([Tol l can be written as N = J n[fjP — Vx 2 , h°]dx, with an 
analogous expression for the density difference. Introducing 
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FIG. 3: (color online) Density and polarization profiles in the asym- 
metric trap for different particle numbers iV = 40, 80, 160 at a fixed 
effective density p e g — N\/V at U = — 8t (we choose V such 
that at each N, p cB = N VV = 20^/a002 =const). (a) p = 0, 
(b) p — 0.5. For the largest N and p = 1/2, we use chains with 
L = 300 sites. 



the new coordinate y = Wx, we get NW = J n\po — 
y, h°], showing that the constants fjP, h° (and through them, 
the density profiles) are completely determined by three pa- 
rameters: the effective density p e g = N^/V, the spin polar- 
ization p = (Nf — N±)/N, and the interaction strength U. 

In this article we mainly consider configurations where the 
core of the cloud is partially polarized, corresponding to tra- 
jectories starting from the PP phase. Then, LDA predicts the 
existence of two shell structures: a partially polarized (PP) 
core with either (i) fully polarized (FPi) or (ii) fully paired 
(ED) wings (see the vertical lines in Fig. Q] for an illustra- 
tion). The dotted line in Fig. Q] corresponds to the special 
case h = y/U 2 + 16/2 — 2, where the PP phase extends to 
the edge of the cloud, in analogy with the continuum case 
EI!!. This yields a critical value of the spin polarization 
p = p c , separating the regime with fully paired (p < p c ) from 
the one with fully polarized wings (p > p c ). At larger densi- 
ties, more shell structures are possible, e.g., a band insulator 
of the majority spins in the center with two surrounding shells 



C. DMRG simulations 

The density profiles of the gas in the trap can be calculated 
to a high accuracy using the DMRG algorithm. We use a stan- 
dard implementation for ground state calculations l44ll and 
up to 600 states. Occasionally, it was necessary to perform 
an unusually large number of sweeps (sometimes more than 
20) to obtain converged results. For a Hamiltonian with site- 
dependent parameters, it is convenient to grow the system in a 
linear fashion, say from the left to the right, during the warm- 
up part of the DMRG algorithm. The largest system size used 
in this work is L = 300 sites for p = 0.5, U = —8t and 
N = 160 particles. In that case, typical discarded weights in 
the center of the system are smaller than 10~ 10 . In the two 
following subsection we provide details on the set-up used in 
the DMRG runs, such as the choice for the trapping potential 
and a finite-size scaling analysis. 
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1. Asymmetric vs symmetric traps 

We here consider two types of a harmonic confining poten- 
tial parameterized by a constant V, namely one that is sym- 
metrically placed in the optical lattice [see Eq. (0] and one 
that traps particles in the left part of the chain: 



trap 



x 2 rii . 



(11) 



We will next show that for the purpose of studying the behav- 
ior in the wings, the symmetric Eq. (0 and the asymmetric 
set-up Eq. (fTTT > yield the same quantitative behavior. Note 
that in the figures, we shall display the results for the symmet- 
ric trap shifted by L/2 lattice sites to the left with respect to 
Eq. (0. The advantage of using the asymmetric trap is that 
it allows us to reach much lower densities and larger particle 
numbers without resorting to large system sizes, and at less 
computational costs. Note that for the same trapping ampli- 
tude V, the particle numbers in the symmetric trap Eq. (0 are 
twice as large as the ones in the corresponding asymmetric 
set-up. 

Figure0shows the comparison between the symmetric and 
the asymmetric trap for the parameters of Ref. QT8Q (U = —8t, 
V = 0.0022;, N = 40 in the case of the symmetric trap) 
for the balanced case [p = 0, panel (a)] and the imbalanced 
case [p = 1/2, panel (b)]. The open boundary that is present 
when we use the asymmetric set-up clearly induces strong ad- 
ditional oscillations in the density profile, both at p = and 
p = 0.5. The main point is, however, that away from the left 
boundary and towards the edge of the cloud, the density and 
polarization profiles computed from either the symmetric or 
the asymmetric trap coincide. In the case of p = 0.5, the sys- 
tem enters the fully polarized wings at the same distance from 
the trap center, independently of the particle number. This jus- 
tifies the use of the asymmetric trap in our analysis of the shell 
structure, which we shall use in the remainder of this work. 



2. Scaling analysis in the particle number 

As a next step we study the finite-size scaling of the den- 
sity and polarization profiles. More precisely, we increase the 
particle number but keep the effective density p ff = N\/V 
fixed ll48l l49ll (the effective density is given for the asymmet- 
ric trap unless stated otherwise). Site-dependent quantities 
should then be plotted versus where £ = (V/t)^ 1 / 2 is 
the characteristic length scale for a given effective density. 
The results of such an analysis are presented in Fig.0for the 
parameters of Fig. [2] with particle numbers N = 40, 80, 160 
in the asymmetric trap. As expected, all curves computed for 
different particle numbers fall on top of each other. We fur- 
ther see that the amplitude of the boundary-induced oscilla- 
tions becomes much weaker the larger the polarization or the 
larger N is and therefore, the presence of the oscillations does 
not pose a problem for our analysis. 
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FIG. 4: (color online) Shell structure in an asymmetric trap for (a) 
U = -8t, N = 160, (b) U = -8t, N = 80 and (c), (d) U = -At, 
N = 160. In all cases, V = 3.125 10~ 5 1 Thick solid lines: density 
profile (m), dashed lines: polarization profile p;, thin solid lines 
in (c), (d): double occupancy 2{di). Panels (a), (c), and (d) show 
DMRG results only while (b) includes LDA results (thin lines) for 
comparison. The inset in (b) shows the density difference and the 
insets in (c), (d) give an enlarged view of the interface regions. 



III. RESULTS 



A. Fully paired wings 



In this section we will focus on N = 
particles, with a trapping potential of V 



160 and N = 
: 3.125 
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change from fully paired (ED, equal densities) wings to fully 
polarized (FP\) wings is expected to occur at small polariza- 



P 



< 



0.2. 



Figure Ua) 
the original 



tions, and we hence concentrate on 
contains the results for p P ff = 20 VO-002, i.e. 
density studied in Ref. 1181 . but with a much larger number 
of particles. While at p = 0.1, we still clearly observe FP 
wings, at p = 0.05, the local polarization practically vanishes 
at the edge of the cloud. Therefore, for these parameters, the 
critical polarization for the change in the shell structure must 
be 0.05 < p c < 0.1, and further analysis shows that it is 
p c Rj 0.07. Therefore, if there are too few particles in the 
trap, say N = 40, then the fully paired wings are very diffi- 
cult to observe since such a p c would correspond to one of the 
smallest possible polarizations. Moreover, effects due to a fi- 
nite particle number become more relevant for small N and/or 
tighter traps, and LDA does not necessarily quantitatively ap- 
ply in this regime (see also Ref. lf22tl ). In fact, numerically ex- 
act DMRG calculations do not yield evidence for fully paired 
wings at, for instance, N = 40 and V = 0.005i lfl8tl , 

By going to half as many particles yet while keeping the 
trapping amplitude constant at V = 3.125 10 -5 i and thus de- 
creasing the effective density, the fully paired wings become 
better visible, as we demonstrate in Fig. |4fb) for p = 0.05. 
There, the local polarization vanishes at a distance from the 
center of the trap close to x/£ ~ 0.5, whereas the cloud ex- 
tends up to x/£ ~ 0.6. Clearly, the volume of the fully paired 
wings increases going to a lower effective density. 

Another way of influencing the volume of the fully paired 
wings is by going to smaller values of the interactions \U\. 
The respective results for U = — At and the same density and 
particle number as in Fig. |4ja) are shown in Figs.|4|c) and (d) 
for p = 0.05 and p — 0.1, respectively. In this example, we 
find evidence for the fully paired wings at p = 0.05, while at 
p = 0.1, we see fully polarized wings. The fully paired wings 
show up as first, a region of practically vanishing local polar- 
ization and second, the double occupancy (di) = (n^n^) 
[thin solid lines in Figs. [4{c) and (d)] remains finite there. 
By contrast, whenever we enter the fully polarized wings, the 
double occupancy goes to zero. Summarizing, so far, we have 
illustrated that the fully paired wings indeed show up in ID 
optical lattices at small polarizations in shallow traps, similar 
to the case of the associated continuum model iflq, ITtL l22ll . 
Decreasing the effective density stabilizes this shell structure 
in the sense that, the relative volume of the fully paired wings 
increases [compare Figs. Ufa) and (b)]. 

Next, we compare the density profiles calculated by DMRG 
and by LDA, where the latter makes use of the exact phase 
diagram for the ID Hubbard model (Fig. [1] see l42tl ). This 
is shown in Fig. |4jb) for U = St and a polarization of 
p = 0.05 at which the wings are fully paired. The agree- 
ment between LDA and DMRG is very convincing which 
is expected in the case of a shallow trap as considered here 
(V = 3.125 10~ 5 i). Most importantly for us here, both meth- 
ods show the fully paired wings and the cloud's spatial ex- 
tension for both the density and density difference agree quite 
well. From the figure, some known shortcomings of LDA are 
obvious: first, LDA fails to reproduce the correct decrease of 
the density profile at the edge of the cloud. This implies that 
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FIG. 5: (color online) Critical polarization p c for the change in the 
shell structure from fully paired wings at p < p c to fully polarized 
wings at p > p c , vs. interaction \U/t\ for two different effective 
densities p eS = N W (V = 3.125 l(T 5 i) with (a) N = 160 and 
(b) iV = 80. The DMRG (squares) results were extracted from runs 
with an asymmetric trap. LDA results are displayed with circles. 
The non monotonic behavior of the critical polarization is a unique 
feature of the lattice, not found for the continuum model. 



within LDA it is possible to simply define the radius of the 
cloud by its extension, whereas for a finite particle number 
this quantity fluctuates. Second, the incommensurate oscilla- 
tions in the density induced by the lattice are not captured by 
LDA either 115011 . While these discrepancies are not crucially 
relevant to the analysis carried out here, there is another fea- 
ture that is absent in the LDA, namely the 2Qfflo oscillation 
in the density difference. This is obvious since LDA assumes 
that the gas is locally homogeneous whereas in DMRG cal- 
culations and in experiments, one deals with a finite particle 
number. 

The 2Qfflo modulation is due to the fact that major- 
ity fermions reside in the nodes of the quasi-condensate 
1 51], which has previously been demonstrated with DMRG 
1 18 , [HI H^[ . In other words, the incommensurate ID FFLO 
state is accompanied with a spin-density wave (see Ref. lf52tl ). 
as is illustrated in the inset of Fig.Ub). This feature, which is 
visible in a trap, where the quasi-condensate is pinned, could 
in principle be observed in experiments with trapped two- 
component Fermi gases J26ll5lll . This highlights the impor- 
tance of using exact numerical techniques such as DMRG and 
QMC that can take the trap into account exactly to comple- 
ment analytical methods. Nevertheless, LDA combined with 
the exact solution of the phase diagram jl^[l7ll has been very 
successful in predicting the shell structures in one dimension, 
which for the continuum case has been observed experimen- 
tally m. 



B. Critical polarization 

We now present our theoretical predictions for the critical 
polarization p c associated to the change in the shell structure 
of the density profiles, from fully paired (p < p c ) to fully 
polarized (p > p c ) wings. For a shallow trap, this quantity 
depends only on the values of the interaction U and the effec- 
tive density p e a = Ny/V. 

In Fig. [5] we plot the critical polarization as a function 
of \U\ comparing DMRG (squares) with LDA (circles), for 
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V = 3.125 10 _5 i and two different particle numbers N = 
160 and N = 80. From the DMRG simulations, we esti- 
mate p c as the point at which, by visual inspection, we can 
clearly resolve regions in the wings with (sj) = Due to 
the finite particle number, these results carry uncertainties of 
Sp c w ±0.0125, 0.03 for N = 80, 160, respectively. 

The agreement between the DMRG results and LDA is ob- 
viously very good. For \ U\ < 5t, the critical polarization is an 
monotonically increasing function of the interaction, in anal- 
ogy with the continuum model lfl6l[T7ll . However, for stronger 
interactions the critical polarization in a lattice reaches a max- 
imum value and then decreases at large \U\, in contrast with 
the continuum case iFIrjfTvTJ . 

The physical explanation is that in this regime the kinetic 
energy of the pairs is suppressed, since their effective mass 
becomes increasingly large (~ \U\/t 2 ) compared to that of 
single fermions (1/t). As a consequence, the energy needed 
to break a pair is given by \U\, independently of the parti- 
cle density. This means that the boundary between the fully 
paired (ED) and the partially polarized phase in Fig.[TJ(the line 
connecting the points A,B, and C counterclockwise) reduces 
to the dotted vertical line for U — > — oo, implying the absence 
of fully paired wings in a trap. 

Figure [5] further suggests that for the chosen parameters, 
the maximum critical polarizations are p™ ax w 0.08 and 
p™ x rj 0.1, respectively. These values are smaller than the 
largest critical polarization p" lax w 0.2 found in the contin- 
uum model. 

In Fig. [6] we plot LDA results for the critical polarization 
versus the effective density for different values of the inter- 



action. We see that the polarization increases by decreasing 
the effective density and saturates at a maximum value of 
p™ ax (C/). Then, if one sends U to zero as well, this maxi- 
mum will approach p™ ax ([/) — > 0.2 as is illustrated in the 
inset of Fig. [6] consistent with the findings reported for the 
continuum limit flrjl [l7tl . 



IV. SUMMARY AND DISCUSSION 

In summary, we investigated the density profiles of trapped 
spin-imbalanced Fermi gases in a ID optical lattice. We 
demonstrated that, when performed on sufficiently large sys- 
tems, numerically exact DMRG simulations are fully consis- 
tent with LDA results. In particular, we reported numerical 
evidence for the fully paired wings that were predicted to ex- 
ist at small polarizations lfl6[[l7ll . However, differently from 
the continuum case, we found that the upper critical polariza- 
tion p c for the appearance of the fully paired wings is a non 
monotonic function of the interaction strength U at a fixed 
filling. This behavior is a direct consequence of the increased 
inertia of the pairs in the strong coupling regime (|f | 3> At). 
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